InĀ [17]:
set.seed(999)
options(scipen = 9)
options(warn = -1) 
source("./environment/libraries.R")
knitr::opts_chunk$set(fig.height = 12, fig.width = 18, fig.dpi = 300)
knitr::opts_chunk$set(warning = FALSE)
All packages loaded successfully
InĀ [18]:
name <- "Kenya_E1"
dataset <- read.csv(paste0("./test/", name, "_processed.csv"))
dataset_c <- data.frame(dataset[7:ncol(dataset)], row.names = dataset$Serial)
proj_coord <- data.frame("Easting" = dataset$Easting, "Northing" = dataset$Northing)
xy_coord <- data.frame("Longitude" = dataset$Longitude, "Latitude" = dataset$Latitude)

dataset_c_closed <- cbind(dataset_c, "Res" = 100 - rowSums(dataset_c)) 
head(dataset_c_closed)

dataset_spdf <- cbind(proj_coord, dataset_c_closed)
rownames(dataset_spdf) <- rownames(dataset_c)

structures_geojson_path <- file.path("./data", paste0(name, "_topo_lines.geojson"))
structures <- st_read(structures_geojson_path, quiet = TRUE)
structures_utm <- st_transform(structures, crs = 32637) # Transform to UTM Zone 35S (WGS84)


output_dirs <- "./test/ck/"
if (!dir.exists(output_dirs)) {
  dir.create(output_dirs, recursive = TRUE)
}
A data.frame: 6 Ɨ 17
PCaTiMnFeNiZnRbSrYZrBaMgAlSiKRes
<dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl>
10.27313.47950.33090.06032.72480.00210.00330.00350.08080.00170.0123000000.09331.2275275.45467821.960071.97141662.32071
20.38353.95870.28280.05662.40540.00170.00310.00370.08390.00160.0041000000.09551.2341004.31907318.230192.22713966.70890
30.41584.07200.22650.06172.16340.00110.00330.00340.09210.00120.0045476330.07551.3168314.05277517.194472.48983467.82554
40.37043.95890.29020.06012.63740.00220.00350.00300.08370.00210.0102000000.08801.3842474.66847818.721382.42592465.29028
50.29813.39790.29840.07072.70030.00260.00300.00350.07970.00170.0152000000.10051.2800484.96270820.899372.48651463.39976
60.43734.00800.18080.05371.92700.00100.00290.00300.08600.00120.0035276360.05831.3045303.23393613.664892.97314372.06078
InĀ [19]:
source("./utils/functions/create_quick_map.R")
options(repr.plot.width = 10, repr.plot.height = 10)

dbscan_result <- dbscan(dataset_spdf[, c("Easting", "Northing")], eps = 10, minPts = 5)
Area <- factor(dbscan_result$cluster)
dataset$Area <- Area
create_quick_map(dataset, structures, group_data = "Area")

dataset_spdf$Area <- Area
dataset_spdf <- dataset_spdf[dataset_spdf$Area == 2, ]
No description has been provided for this image
InĀ [20]:
create_spatial_objects <- function(data_spdf, transformation = "none") {
  coords <- data_spdf[, c("Easting", "Northing")]
  crs_string <- "+proj=utm +zone=37 +north +datum=WGS84"
  
  if (transformation == "ilr") {
    spatial_data <- as.data.frame(ilr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
  } else if (transformation == "clr") {
    spatial_data <- as.data.frame(clr(data_spdf[, -c(1, 2, ncol(data_spdf))]))
  } else {
    spatial_data <- data_spdf[, -c(1, 2, ncol(data_spdf))]
  }
  
  return(SpatialPointsDataFrame(coords = coords, data = spatial_data, 
                               proj4string = CRS(crs_string)))
}

spdf_comp <- create_spatial_objects(dataset_spdf, "none")
spdf_ilr <- create_spatial_objects(dataset_spdf, "ilr")
spdf_clr <- create_spatial_objects(dataset_spdf, "clr")

source("./utils/functions/lag_distance_from_spdf.R")
source("./utils/functions/site_diagonal_from_spdf.R")
lag_dist <- lag_distance_from_spdf(spdf_ilr)
site_diag <- site_diagonal_from_spdf(spdf_ilr)
InĀ [21]:
source("./utils/functions/geostatistical_analysis.R")
source("./utils/functions/validate_geostatistical_predictions.R")
InĀ [22]:
results_ilr_omni <- geostatistical_analysis(spdf_ilr, "ILR_Omnidirectional", lag_dist, site_diag, directional = FALSE)
Selected models:
          psill    range kappa   n
Nug 0.005636354 0.000000     0 136
Sph 0.013267064 4.879795     0  69
Gau 0.006745542 3.073336     0  38
Exp 0.018271788 1.914671     0  29
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-omnidirectional/maps/ 
No description has been provided for this image
InĀ [23]:
validate_predictions(spdf_ilr, results_ilr_omni$ck, results_ilr_omni$g)
               V1            V2            V3            V4            V5
ME   0.0001022348 -0.0062020027 -0.0013457263 -0.0008333475  0.0007015358
MSE   0.008724575   0.072485737   0.018183152   0.013407582   0.037057207
               V6            V7            V8            V9           V10
ME   0.0015601790  0.0031828449  0.0014674274 -0.0015878812 -0.0082124051
MSE   0.008616293   0.021176479   0.027666382   0.018340586   0.103614886
              V11           V12           V13           V14           V15
ME   0.0014876697  0.0015760092  0.0003796410  0.0011482589  0.0017138209
MSE   0.023532552   0.012736559   0.021056677   0.027255291   0.011694586
              V16
ME   0.0007488469
MSE   0.010896231
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
$mv_results
A data.frame: 1 Ɨ 3
AccuracyPrecisionGoodness
<dbl><dbl><dbl>
0.850.640.82
$var_results
A data.frame: 16 Ɨ 4
VariableAccuracyPrecisionGoodness
<chr><dbl><dbl><dbl>
V1 0.950.860.93
V2 0.850.780.89
V3 0.900.870.93
V4 0.950.790.89
V5 0.950.930.97
V6 0.850.950.97
V7 0.900.900.95
V8 0.950.830.92
V9 0.950.870.93
V100.900.900.95
V110.900.860.93
V120.950.840.92
V130.950.780.89
V140.950.800.90
V150.950.830.92
V160.950.910.95
No description has been provided for this image
InĀ [24]:
results_ilr_dir <- geostatistical_analysis(spdf_ilr, "ILR_Directional", lag_dist, site_diag, directional = FALSE)
Selected models:
          psill    range kappa   n
Nug 0.005636354 0.000000     0 136
Sph 0.013267064 4.879795     0  69
Gau 0.006745542 3.073336     0  38
Exp 0.018271788 1.914671     0  29
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
Saved 17 maps to: ./test/ck/ilr-directional/maps/ 
No description has been provided for this image
InĀ [25]:
validate_predictions(spdf_ilr, results_ilr_dir$ck, results_ilr_dir$g)
               V1            V2            V3            V4            V5
ME  -0.0005531453  0.0048508733  0.0024681402  0.0026194095  0.0043327897
MSE   0.008576134   0.078193724   0.018882226   0.013745652   0.038442514
               V6            V7            V8            V9           V10
ME   0.0016470742  0.0017199244 -0.0029722507  0.0016671107  0.0004749151
MSE   0.008548680   0.022354672   0.028832000   0.017899969   0.101678939
              V11           V12           V13           V14           V15
ME   0.0005650498  0.0016746456  0.0042720624  0.0037208622  0.0015264638
MSE   0.024453100   0.013906207   0.019529877   0.026214458   0.012438465
              V16
ME  -0.0006585043
MSE   0.011510467
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
$mv_results
A data.frame: 1 Ɨ 3
AccuracyPrecisionGoodness
<dbl><dbl><dbl>
0.850.640.82
$var_results
A data.frame: 16 Ɨ 4
VariableAccuracyPrecisionGoodness
<chr><dbl><dbl><dbl>
V1 0.950.860.93
V2 0.850.780.89
V3 0.900.870.93
V4 0.950.790.89
V5 0.950.930.97
V6 0.850.950.97
V7 0.900.900.95
V8 0.950.830.92
V9 0.950.870.93
V100.900.900.95
V110.900.860.93
V120.950.840.92
V130.950.780.89
V140.950.800.90
V150.950.830.92
V160.950.910.95
No description has been provided for this image
InĀ [26]:
par(bg = "white")
options(repr.plot.width = 15, repr.plot.height = 12)

# Compute MAF analysis
g_temp <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_omni <- variogram(g_temp, width = lag_dist / 2, cutoff = site_diag / 3, cross = TRUE)
maf_result <- Maf(spdf_ilr@data, vg = v_omni)

source("./utils/functions/quick_pca_screeplot.R")
quick_pca_screeplot(maf_result)
maf_dataset <- maf_result$scores[, 1:9]
spdf_maf <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                  data = as.data.frame(maf_dataset),
                                  proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))

results_maf_omni <- geostatistical_analysis(spdf_maf, "MAF_Omnidirectional", lag_dist, site_diag,
                                           directional = FALSE, maf_result = maf_result, 
                                           resolution = 2)
Selected models:
        psill    range kappa  n
Nug 0.4284907 0.000000     0 45
Exp 0.4890494 1.198881     0 24
Sph 0.7692933 3.407950     0 21
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
No description has been provided for this image
Saved 17 maps to: ./test/ck/maf-omnidirectional/maps/ 
No description has been provided for this image
InĀ [27]:
validate_predictions(spdf_maf, results_maf_omni$ck, results_maf_omni$g)
            maf1         maf2         maf3         maf4         maf5
ME   0.007983487  0.009709664  0.015301634  0.012281329  0.007186930
MSE    0.9761022    1.0474878    0.9136258    0.9094230    0.8691944
            maf6         maf7         maf8         maf9
ME  -0.008298787  0.021281565 -0.007840009 -0.004253346
MSE    0.8925875    0.6713064    0.7971482    0.9827595
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
$mv_results
A data.frame: 1 Ɨ 3
AccuracyPrecisionGoodness
<dbl><dbl><dbl>
0.950.540.77
$var_results
A data.frame: 9 Ɨ 4
VariableAccuracyPrecisionGoodness
<chr><dbl><dbl><dbl>
maf10.900.930.97
maf20.950.790.89
maf30.950.910.95
maf40.950.830.91
maf50.650.940.96
maf60.950.890.94
maf70.900.810.91
maf80.950.580.79
maf90.900.920.96
No description has been provided for this image
InĀ [28]:
g_temp2 <- create_gstat_from_spdf(spdf_ilr, method = "ordinary")
v_dir <- variogram(g_temp2, width = lag_dist / 2, cutoff = site_diag / 3,
                  alpha = c(0, 45, 90, 135), tol.hor = 22.5, cross = FALSE)

maf_result2 <- Maf(spdf_ilr@data, vg = v_dir)
quick_pca_screeplot(maf_result2)
maf_dataset2 <- maf_result2$scores[, 1:10]
spdf_maf2 <- SpatialPointsDataFrame(coords = dataset_spdf[, c("Easting", "Northing")], 
                                   data = as.data.frame(maf_dataset2),
                                   proj4string = CRS("+proj=utm +zone=37 +north +datum=WGS84"))

results_maf_dir <- geostatistical_analysis(spdf_maf2, "MAF_Directional", lag_dist, site_diag,
                                          directional = TRUE, maf_result = maf_result2, 
                                          resolution = 2)
Directional approach not fully implemented yet - using omnidirectional as fallback
Selected models:
        psill    range kappa  n
Nug 0.5618877 0.000000     0 55
Sph 0.4160659 3.029266     0 34
Gau 0.9858814 1.388666     0 13
Exp 0.5556532 1.361615     0  8
Linear Model of Coregionalization found. Good.
[using ordinary cokriging]
No description has been provided for this image
Saved 17 maps to: ./test/ck/maf-directional/maps/ 
No description has been provided for this image
InĀ [29]:
validate_predictions(spdf_maf2, results_maf_dir$ck, results_maf_dir$g)
            maf1         maf2         maf3         maf4         maf5
ME   0.008505926  0.010039250  0.004414058 -0.007083254 -0.010341529
MSE    0.9789648    0.9862980    1.0416295    0.9492662    0.8339911
            maf6         maf7         maf8         maf9        maf10
ME  -0.010617249  0.021731741  0.006564275 -0.001078079  0.014045939
MSE    0.9829563    0.9810363    0.6349502    0.9448917    0.7574230
No description has been provided for this image
No description has been provided for this image
No description has been provided for this image
$mv_results
A data.frame: 1 Ɨ 3
AccuracyPrecisionGoodness
<dbl><dbl><dbl>
0.950.490.74
$var_results
A data.frame: 10 Ɨ 4
VariableAccuracyPrecisionGoodness
<chr><dbl><dbl><dbl>
maf1 0.950.930.96
maf2 0.800.960.98
maf3 0.950.830.92
maf4 0.950.900.95
maf5 0.950.760.88
maf6 0.900.890.94
maf7 0.950.920.96
maf8 0.950.770.88
maf9 0.950.650.83
maf100.950.810.91
No description has been provided for this image
InĀ [30]:
source("./utils/functions/create_ck_maplist_from_list.R")

ck_files <- list.files("./test/ck/", full.names = TRUE, pattern = paste0(name, ".*\\.rds$"), recursive = TRUE)
ck_list <- lapply(ck_files, readRDS)

if (length(ck_list) > 0) {
  cat("Available cokriging results:", length(ck_list), "\n")
  cat("Files:", basename(ck_files), "\n")
  
  if (length(ck_list) > 1) {
    cat("Multiple methods completed. Individual results saved in method subdirectories\n")
  }
} else {
  cat("No cokriging results found\n")
}
cat("Geostatistical workflow completed successfully!\n")
Available cokriging results: 4 
Files: Kenya_E1_CK2.rds Kenya_E1_CK2.rds Kenya_E1_CK2.rds Kenya_E1_CK2.rds 
Multiple methods completed. Individual results saved in method subdirectories
Geostatistical workflow completed successfully!